#include<Algorithm/fast_compute_maybe_intersecting_areas.h>

#include<tool/libcell_tools_octree.h>
#include<tool/libcell_octree_iterator.h>
#include<Mesh/libcell_iterator.h>
#define DEFAULT_THRESHOLD  1e-6 
#define quote lib_cell_quote



static inline void get_cube_from_cell(template_c*c,double* re)
{
    //double* re=(double*)malloc(sizeof(double)*6);
    memset(re,0,sizeof(double)*6);
    re[0]=c->vertices[0]->point[0];re[1]=c->vertices[0]->point[1];re[2]=c->vertices[0]->point[2];
    re[3]=c->vertices[0]->point[0];re[4]=c->vertices[0]->point[1];re[5]=c->vertices[0]->point[2];
    for(int i=1;i<c->vertices_size;i++)
    {
        for(int j=0;j<3;j++)
        {
            c->vertices[i]->point[j]<re[j]?(re[j]=c->vertices[i]->point[j]):(c->vertices[i]->point[j]>re[j+3]?(re[j+3]=c->vertices[i]->point[j]):i=i);
        }
    }
    //return re; 
}



static inline int is_intersected_two_cube(double *min1,double* max1,double* min2,double* max2)
{
    if(min2[0]>max1[0]+DEFAULT_THRESHOLD||min2[1]>max1[1]+DEFAULT_THRESHOLD||min2[2]>max1[2]+DEFAULT_THRESHOLD)
    {
        return 0;
    }
    if(max2[0]<min1[0]+DEFAULT_THRESHOLD||max2[1]<min1[1]+DEFAULT_THRESHOLD||max2[2]<min1[2]+DEFAULT_THRESHOLD) 
    {
        return 0;
    }
    return 1;
}

static Node* get_vertices_from_node_of_cells(Node* n)
{
    Node* re=NULL;

    Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        for(int i=0;i<c->vertices_size;i++)
        {
            if(tree1->find(tree1,c->vertices[i]->id)==NULL)
            {
                tree1->insert(tree1,c->vertices[i]->id,c->vertices[i]);
                re=node_overlying(re,c->vertices[i]);
            }
        }
    } 
    int_rb_tree_free(tree1);
    return re; 
}

static double * get_bounding_box_from_nodes_of_cells(Node* n)
{
    if(n==NULL)
    {
        return NULL;
    }
    double *re=(double*)malloc(sizeof(double)*6);
    Node* nit=n;
    template_c* c=(template_c*)(nit->value);
    re[0]=c->vertices[0]->point[0]; re[1]=c->vertices[0]->point[1];  re[2]=c->vertices[0]->point[2];
    re[3]=c->vertices[0]->point[0]; re[4]=c->vertices[0]->point[1];  re[5]=c->vertices[0]->point[2];
    for(;nit!=NULL;nit=(Node*)(nit->Next))
    {
        c=(template_c*)(nit->value);
        for(int j=0;j<c->vertices_size;j++)
        {
            template_v* v=c->vertices[j];
            for(int i=0;i<3;i++)
            {
                v->point[i]>re[3+i]?(re[3+i]=v->point[i]):(v->point[i]<re[i]?(re[i]=v->point[i]):(i=i));
            } 
        }
    }
    if(re[3]-re[0]<=1e-3){re[3]+=0.001;re[0]-=0.001;}
    if(re[4]-re[1]<=1e-3){re[4]+=0.001;re[1]-=0.001;}
    if(re[5]-re[2]<=1e-3){re[5]+=0.001;re[2]-=0.001;}

    return re;
}



static Node* simplify_node_of_two_nodes(Node* n1,Node *n2,int* flag)
{
    *flag=0;
    if(n1==NULL||n2==NULL)
    {
        return NULL;
    }
    Node* re=NULL;
    double* cube1=get_bounding_box_from_nodes_of_cells(n1);
    double cube2[6]={0};
    for(Node* nit=n2;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        get_cube_from_cell(c,cube2);
        if(is_intersected_two_cube(cube1,&(cube1[3]),cube2,&(cube2[3])))
        {re=node_overlying(re,c);}
        else
        {*flag=1;}
        
    }
    free(cube1);
    return re;
}


static double * get_bounding_box_from_nodes_of_vertices(Node*n )
{
    if(n==NULL)
    {
        return NULL;
    }
    double *re=(double*)malloc(sizeof(double)*6);
    Node* nit=n;
    template_v* v=(template_v*)(nit->value);
    re[0]=v->point[0]; re[1]=v->point[1];  re[2]=v->point[2];
    re[3]=v->point[0]; re[4]=v->point[1];  re[5]=v->point[2];

    for(;nit!=NULL;nit=(Node*)(nit->Next))
    {
        v=(template_v*)(nit->value);
        for(int i=0;i<3;i++)
        {
            v->point[i]>re[3+i]?(re[3+i]=v->point[i]):(v->point[i]<re[i]?(re[i]=v->point[i]):(i=i));
        } 
    }
    if(re[3]-re[0]<=1e-4){re[3]+=0.0001;re[0]-=0.0001;}
    if(re[4]-re[1]<=1e-4){re[4]+=0.0001;re[1]-=0.0001;}
    if(re[5]-re[2]<=1e-4){re[5]+=0.0001;re[2]-=0.0001;}

    return re;
}

typedef struct  Temp_Info{
    Node* cells;
    double * bounding_box;
}Temp_Info;
static inline  void temp_info_init(Temp_Info* ti)
{
    ti->cells=NULL;
    ti->bounding_box=NULL;
}
static inline void temp_info_free(Temp_Info*ti)
{
    if(ti==NULL)
    {
        return;
    }
    free_node(ti->cells);
    SAFE_FREE(ti->bounding_box);
    free(ti);
}


static Temp_Info* my_get_temp_info(Node* n)
{
    Temp_Info* ti=(Temp_Info*)malloc(sizeof(Temp_Info));
    temp_info_init(ti);
    Int_RB_Tree*tree2=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree2);
    Node* nit=n; 
    template_v* v=(template_v*)(nit->value);
    
    for(;nit!=NULL;nit=(Node*)(nit->Next))
    {
        v=(template_v*)(nit->value);
        
        // for(Node* vcit=v->cells;vcit!=NULL;vcit=(Node*)(vcit->Next))
        for(int i=0;i<v->cellsn.size;i++)
        {
            template_c* c=(template_c*)(v->cellsn.extra_data[i]);
            if(tree2->find(tree2,c->id)==NULL)
            {
                tree2->insert(tree2,c->id,c);
                ti->cells=node_overlying(ti->cells,c);
            }  
        } 
    }  
    int_rb_tree_free(tree2);
    ti->bounding_box=get_bounding_box_from_nodes_of_cells(ti->cells);
    return ti;
}



void simplify_maybe_intersecting_areas(Maybe_Intersecting_Areas* mia)
{
    Node* n1=mia->node1,*n2=mia->node2;
    int flag=0,flag1=0;
    while(1)
    {
        flag=0;
        Node* temp_n1=simplify_node_of_two_nodes(n1,n2,&flag1);
        if(flag1==0)
        {
            flag++;
        }
        free_node(n2);
        n2=n1;n1=temp_n1;temp_n1=NULL;
        temp_n1=simplify_node_of_two_nodes(n1,n2,&flag1);
        if(flag1==0)
        {
            flag++;
        } 
        free_node(n2);
        n2=n1, n1= temp_n1;temp_n1=NULL;
        if(flag>=2)
        {
            break;
        }
    }  
    mia->node1=n1;mia->node2=n2; 
}

Node* get_maybe_intersecting_areas_from_nodes(Node*n1,Node*n2)
{
    OC_Node* ocn1=(OC_Node*)malloc(sizeof(OC_Node)),*ocn2=(OC_Node*)malloc(sizeof(OC_Node)); 
    Node* nn1=get_vertices_from_node_of_cells(n1),*nn2=get_vertices_from_node_of_cells(n2);

    oc_node_init_from_node(ocn1,nn1);oc_node_init_from_node(ocn2,nn2);
    int len1=node_size(nn1),len2=node_size(nn2);
    int ceng1=(len1<4000?2:(len1<9000?3:(len1<20000?4:(len1<100000?5:6) )));
    int ceng2=(len2<4000?2:(len2<9000?3:(len2<20000?4:(len1<100000?5:6) )));
    printf("begin ocdivide\n");
    for(int i=0;i<6;i++)
    {
        if(i<ceng1)
        {oc_node_divide_all_leaves(ocn1);}
        if(i<ceng2)
        {oc_node_divide_all_leaves(ocn2);}
    }
    printf("end ocdivide\n"); 
    Node* tin1=NULL,*tin2=NULL; 
    for(auto ocit=octree_begin(ocn1);ocit.it!=NULL;ocit++)
    {
        if(ocit.it->value!=NULL)
        {
            tin1=node_overlying(tin1,my_get_temp_info(ocit.it->value));
            free_node(ocit.it->value);ocit.it->value=NULL;
        }
    }
    for(auto ocit=octree_begin(ocn2);ocit.it!=NULL;ocit++)
    {
        if(ocit.it->value!=NULL)
        {
            tin2=node_overlying(tin2,my_get_temp_info(ocit.it->value));
            free_node(ocit.it->value);ocit.it->value=NULL;
        }
    }
    Node* re=NULL;
   
    for(Node* nit1=tin1;nit1!=NULL;nit1=(Node*)(nit1->Next))
    {
        Temp_Info* ti1=(Temp_Info*)(nit1->value); 
        Maybe_Intersecting_Areas* mia=(Maybe_Intersecting_Areas*)malloc(sizeof(Maybe_Intersecting_Areas));
        maybe_intersecting_areas_init(mia);
        mia->node1=node_copy(ti1->cells);
        for(Node* nit2=tin2;nit2!=NULL;nit2=(Node*)(nit2->Next))
        {
            Temp_Info* ti2=(Temp_Info*)(nit2->value);
            if(ti2->cells!=NULL&&is_intersected_two_cube(ti1->bounding_box,&(ti1->bounding_box[3]),ti2->bounding_box,&(ti2->bounding_box[3]))==1)
            {
                Node* temp_n=node_copy(ti2->cells);
                Node* temp_nr=node_reverse(temp_n);
                temp_nr->Next=mia->node2;
                if(mia->node2!=NULL)
                {
                    mia->node2->Prev=temp_nr;
                }
                mia->node2=temp_n;
            }
        }

        simplify_maybe_intersecting_areas(mia);
        if(mia->node2==NULL)
        {
             maybe_intersecting_areas_free(mia);
        }  
        else
        {
            re=node_overlying(re,mia);
        }
    } 
    for(Node* nit=tin1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        Temp_Info* ti=(Temp_Info*)(nit->value);
        temp_info_free(ti); 
    } 
    for(Node* nit=tin2;nit!=NULL;nit=(Node*)(nit->Next))
    {
        Temp_Info* ti=(Temp_Info*)(nit->value);
        temp_info_free(ti); 
    } 
    free_node(tin1);free_node(tin2); 
    oc_node_free(ocn1);oc_node_free(ocn2);
    free_node(nn1);free_node(nn2);
    return re;
}

